#!/usr/bin/perl -w
use strict;
use FindBin;
use Net::FTP;
use LWP::Simple;
use HTTP::Date;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin");
use run_cmd;
use run_db_cmd;
use dbh_helper;
my ($cmd);

print STDERR <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    table_cdna_sequence

	Description:
		Maps cDNA sequence against peptide sequence
		Enters both peptide and cDNA sequence into database
		Uses the latter to fill in the pi field in the
			ens_id field
	Populates DB tables:
		sequences
		accessions
		taxon.cdna_sequence
		ens_id

888888888888888888888888888888888888888888888888888888888888888888888

HEADLINE

#-----------------------------------------------------------------------------------------
#	Get DB parameters
#
use db_parameters;
use ortho_species;
use ensembl_parameters;
my @db_param = get_ensembl_taxon_data();
#
#	Get DB parameters
#-----------------------------------------------------------------------------------------

my ($curr_dir) = $FindBin::Bin;
$curr_dir =~ s/^(.*)\/?$/$1/;
my $temp_dir = $dir_pipeline_temp;


my $dbh = connect_to_panda();
use get_psql_time;
my $start_time = get_psql_time($dbh);

#
# get ENST to ENSP mapping
#
print "\tGet ENST to ENSP mapping...\n";
run_psql_cmd("SELECT transcript_id, prot_id from ens_id",
					"$temp_dir/ENST_to_ENSP.map");

unlink ("$temp_dir/aligned.cdna.table.data");

my %taxa_data = get_taxon_data();
#_________________________________________________________________________________________
#	for each taxon, get cDNA and peptide sequence from ENSEMBL
#_________________________________________________________________________________________
for my $taxon(keys %taxa_data)
{
	my $sequence_dir = $taxa_data{$taxon}[SEQ_DIR];

	my $prot_file	= "$sequence_dir/prot.fa";
	my $cdna_file	= "$sequence_dir/cdna.fa";


	print "\n$taxon:\n\n";
	print "\t\tPrepare $taxon cDNA and peptide sequence files for panda:\n";




	#
	# Replace ENST with ENSP
	#
	print "\t\tReplace ENST with ENSP...\n";
	my $cmd = <<"    CMD";
	$dir_lg_share_bin/re_replace_filter_with_lookup.pl
		--from_regex "^>([^ \\\\t]+)"
		--to_regex ">\\\$match \\\$1"
		--lookup_file $temp_dir/ENST_to_ENSP.map
		--keep_non_matching_lines
			$cdna_file	> $temp_dir/cdna.data
    CMD
	#print STDERR $cmd;
	run_cmd($cmd);



	#
	# count the number of stop codons in the peptide file
	#
	$cmd = <<"    CMD"; $cmd =~tr/\t\n/ /d; $cmd =~ s/\s+/ /gm;
	$dir_lg_share_bin/fa_count_stop_codons.pl
		< $prot_file
		$dir_pipeline_errors/table_sequences.$taxon.disrupted_sequences.errors
    CMD
	my $cnt_stops = `$cmd`;
	if ($cnt_stops)
	{
		print "\n\n\t", "!" x 62, "\n\tWarning!\n",
					 "\t\t$cnt_stops\tpeptides with stops found.\n\t", "!" x 62, "\n\n";
	}

	#
	# map protein to cDNA
	#
	print "\tMake sure cDNA sequences and peptides match exactly.\n";
	$cmd = <<"    CMD";
    $dir_lg_share_bin/map_cdna_onto_aa
        --x_for_aa_mismatches   X
		--peptide_sequences		$prot_file
		--cdna_sequences		$temp_dir/cdna.data
		--peptide_output  		$temp_dir/aligned.pep.$taxon.data
		--cdna_output  			$temp_dir/aligned.cdna.$taxon.data
		--table_of_mapped_cdna  $temp_dir/aligned.cdna.$taxon.table.data
		--regex					"[^ \\t]+"
		--err_log				$dir_pipeline_errors/table_sequences.$taxon.cDNA_to_AA_map_errors
		--missing_identifiers	$dir_pipeline_errors/table_sequences.$taxon.cDNA_not_in_prot_errors
		--matching              4
		--verbose
    CMD
	#print STDERR $cmd;
	run_cmd($cmd);

	run_cmd("cat $temp_dir/aligned.cdna.$taxon.table.data >> $temp_dir/aligned.cdna.table.data");


	#
	# panda
	#
	my $db_origin		= $taxa_data{$taxon}[PANDA_DB_ORIGIN];
	print "\t\tEnter peptides into panda...\n";
	$cmd = <<"    CMD";
	$dir_lg_share_bin/panda
			--dbname	$db_name
			--user		$db_user
			--host		$db_host
			--fasta_file	$temp_dir/aligned.pep.$taxon.data
			--name_source "$db_origin"
			--err_log     $dir_pipeline_errors/table_sequences.$taxon.panda_errors
			--verbose
			--taxid $taxa_data{$taxon}[TAX_ID]
    CMD
	$cmd =~tr/\t\n/ /d;
	$cmd =~ s/\s+/ /gm;
	#print STDERR "$cmd\n";
	run_cmd($cmd);

}


	#
	# convert hash to pi
	#
	print "Add cDNA sequences linked to PI in panda...\n";
	$cmd = <<"    CMD";
    $dir_lg_share_bin/add_sequence_linked_data
            --dbname    $db_name
            --user      $db_user
            --host      $db_host
            --err_log   $dir_pipeline_errors/table_sequences.add_cdna_errors
            --verbose
            --table "cdna_sequences"
            --columns "prot_id TEXT, cdna TEXT, cdna_hash OID[]"
            --remove_first
            --input_file $temp_dir/aligned.cdna.table.data
    CMD
	#print STDERR $cmd;
	run_cmd($cmd);

# make sure cdna_sequences unique
$dbh->do("SELECT DISTINCT ON (prot_id) * into temp t_cdna_sequences from cdna_sequences;")
	or die $dbh->errstr."."; 
	
$dbh->do("TRUNCATE cdna_sequences;")
	or die $dbh->errstr."."; 

$dbh->do("INSERT INTO cdna_sequences select * from t_cdna_sequences;")
	or die $dbh->errstr."."; 

$dbh->do("VACUUM ANALYSE cdna_sequences")
	or die $dbh->errstr."."; 


#
# Reindex
#
print "\tReindex cDNA tables...\n";
$dbh->do(<<'PL/SQLCMD')	or die $dbh->errstr.".";
    set SORT_MEM=128000;
	REINDEX INDEX i__cdna_sequences1;
    REINDEX INDEX i__cdna_sequences2;
PL/SQLCMD


#_________________________________________________________________________________________

#	Cleanup and finish

#_________________________________________________________________________________________
unlink glob("$temp_dir/*.data");
unlink glob("$temp_dir/*.map");

#remove zero size (empty error/status) files
(-z $_) && unlink $_ for (glob("$dir_pipeline_errors/*"));
(-z $_) && unlink $_ for (glob("$temp_dir/*"));
#print STDERR "\n";
log_pipeline_stage($dbh, 'NULL', 21, 'table sequences', $start_time);
$dbh->disconnect();
